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Abstract. In a recent publication, we proposed that inflationary perturbation theory 
can be reformulated in terms of a probability transport equation, whose moments 
determine the correlation properties of the primordial curvature perturbation. In this 
paper we generalize this formulation to an arbitrary number of fields. We deduce 
ordinary differential equations for the evolution of the moments of £ on superhorizon 
scales, which can be used to obtain an evolution equation for the dimensionless 
bispectrum, /nl- Our equations are covariant in field space and allow identification 
of the source terms responsible for evolution of /nl- In a model with M scalar 
fields, the number of numerical integrations required to obtain solutions of these 
equations scales like 0(M 3 ). The performance of the moment transport algorithm 
means that numerical calculations with M 3> 1 fields are straightforward. We illustrate 
this performance with a numerical calculation of /nl in Nflation models containing 
M ~ 10 2 fields, finding agreement with existing analytic calculations. We comment 
briefly on extensions of the method beyond the slow-roll approximation, or to calculate 
higher order parameters such as ^nl- 
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1. Introduction 

In a previous article pQ we introduced a new technique ( "moment transport" ) designed 
to extract the statistics of primordial fluctuations generated by inflation. Our technique 
exploits the principle that one never uses a model of inflation to predict the properties 
of any particular universe. Instead, one determines averages over a large ensemble 
of universes. Observable quantities are predicted on the basis that our universe is a 
typical member of the ensemble. Therefore there is no need to study the evolution 
of fluctuations in any specific universe; it is sufficient to compute the evolution of the 
probability density function on the ensemble. Ref. pQ discussed the moment-transport 
method for a restricted range of models. In this paper we refine and extend our method 
to accommodate models containing many self-interacting scalar fields in a variety of 
gauges. We show that the moment-transport technique yields a system of coupled 
ordinary differential equations which can be used for rapid computation of observable 
quantities, such as the dimensionless bispectrum /nl- 

Our method is based on a transport equation. This is a partial differential equation 
which governs the evolution of a probability distribution. In Ref. pQ, we showed how a 
transport equation for a bivariate distribution can be used to derive ordinary differential 
equations which describe the evolving moments of any nearly Gaussian distribution. Our 
method was based on diagonalization of the covariance matrix describing correlations 
between the two variables, a process analogous to Gram-Schmidt orthogonalization. In 
this paper our first objective is to introduce a more elegant method based on matrix 
decomposition. Such decompositions generalize immediately to an arbitrary number of 
scalar fields. Using the resulting equations we are able to evolve the joint probability 
distribution for any number of variables. 

The moment transport equation approach can be thought of as a reformulation of 
cosmological perturbation theory on superhorizon scales with the additional virtue that, 
rather than simply evolving a perturbed quantity, such as (, this approach evolves the 
statistical moments of a perturbed quantity. This is advantageous because these are the 
observationally relevant objects. 

The equations presented in Ref. pQ were valid only in the spatially flat slicing. 
In this gauge, the joint probability distribution for a collection of light scalar fields 
at horizon- crossing is simple to calculate using the methods of quantum field theory 
[2] . Ref. pQ determined the statistics of observational quantities — such as the curvature 
perturbation, (, in the uniform density gauge — by evolving this distribution to a later 
time and making an appropriate gauge transformation. In the present paper, our 
second objective is to generalize the moment-transport method to arbitrary slicings 
of spacetime. Once liberated from the restriction to spatially flat slicings we are free 
to calculate the horizon- crossing distribution of the curvature perturbation itself, and 
propagate this forward in time. Observational quantities may be read directly from this 
distribution at any desired time, with no need for an auxiliary gauge transformation. 

Whichever gauge we pick, our predictions for observational quantities must agree. 
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Nevertheless, working directly in the uniform density gauge has certain advantages. 
First, slices of constant time have an intuitive and unambiguous meaning, being 
associated with hypersurfaces of constant Hubble parameter. Second, it is possible 
to write evolution equations for observable quantities such as /nl- There is no need 
to keep track of intermediate, unphysical quantities with the associated risk that large 
cancellations occur when taking combinations designed to yield physical observables. 
Third, these evolution equations show explicitly which source terms are responsible for 
the growth and decay of the covariance and higher moments of (. The possibility of 
cancellations obscures this interpretation for intermediate variables such as the field 
perturbations. 

In this paper, we discuss two methods which may be used to extract the moment 
hierarchy associated with a general transport equation. In $2] we revisit the technique of 
Gauss-Hermite expansions, employed in Ref. [1], which explicitly invokes a perturbative 
expansion around a Gaussian distribution. Similar expressions were used by Contaldi 
& Magueijo to synthesize microwave background maps with a nongaussian component 
[3], and more recently have been applied to the distribution of collapsed structures; 
among others, see Refs. [U EJ E]. Improving the discussion given in Ref. P, we use a 
technique of matrix decomposition to find equations valid for an arbitrary number of 
fields. This method has several virtues. Most important, in the generalized version to 
be discussed in §2] it is explicitly tensorial. Therefore, combinatorical factors associated 
with multi-field models are built into the formalism and do not need to be addressed 
directly. 

The Gauss-Hermite technique is an example of a cumulant expansion. Given a 
"kernel" distribution with a finite number of nonzero cumulants, up to some order M, 
we use it as a template for a general distribution with perturbatively small cumulants 
of order > M. For the Gauss-Hermite expansion the kernel is a Gaussian with M = 2. 
It has zero third- and higher-order cumulants. In the absence of special reasons to 
the contrary, the Gauss-Hermite method fails when a third- or higher cumulant of the 
probability density function of interest grows to the degree that it is not perturbatively 
small. 

In $3] we give a different perspective on the results of ^2j rederiving them using an 
alternative technique which does not rely on an expansion around an unperturbed kernel. 
Introducing generating functions for the moments and cumulants of the probability 
density function of interest, we demonstrate the surprising and remarkable fact that the 
same transport hierarchy applies for an arbitrary kernel function. The disadvantage 
of this method is a complicated treatment of the multi-field combinatorics. As well 
as a check on the correctness of the formulas derived in $31 the method of generating 
functions clarifies under which circumstances we can expect the transport hierarchy 
to apply. In particular, it allows a more refined discussion of the error involved in 
truncation of the hierarchy. In practical calculations the two methods are essentially 
equivalent because the transport hierarchy must be truncated by assuming that an 
infinite number of cumulants are perturbatively small. 
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In §H we apply our general framework to the fluctuations generated by a model of 
inflation. This depends on the selection of a smoothing scale and a choice of gauge, 
which we discuss in §4.11 The smoothing scale is a measure of the spatial scale on 
which the probability density function measures correlations. In §4.21 we specialize to 
the uniform density gauge and obtain the moment transport equations on this slicing. 

Throughout, we adopt units in which c = h~ = 1 but explicitly retain the 
Planck mass, Mp = (87rG) -1 / 2 , because in some gauges we are obliged to mix 
dimensionless quantities — such as the accumulated e-folds, N, and its perturbation (— 
with dimensionful ones, such as the field values 0,. We label the species of scalar fields 
by indices 

2. Method A: Gauss— Hermite expansions 

In this section, we use the Gauss-Hermite method to derive the moment hierarchy 
associated with a transport equation for an arbitrary number of variables a?j, where 
% — 1, . . . , N. In due course these variables will be asociated with cosmological quantities 
such as field values or curvature, but at present we keep the discussion general. It will 
sometimes be convenient to adopt a matrix notation, in which the X\ are treated as 
components of a vector x. In what follows we use matrix and component notation 
interchangeably. 

2.1. The Gauss-Hermite expansion 

We allow x to have an arbitrary expectation value X(£). Two-point correlations among 
the x are expressed by the covariance matrix which is defined to satisfy 



The covariance matrix and third-order moments are functions of time, as is the position 
of the centroid X. As in Eqs. ([1]) — d2J) , we will frequently suppress explicit time 
dependence to avoid unnecessary clutter. 

Eq. (pQ) makes S a real, symmetric positive-definite matrix. It therefore admits a 
decomposition of the form X = AA T , where A T denotes the matrix transpose of A. A 
candidate decomposition can be obtained by setting A = QX 1 ^ 2 where A = Q T HQ is 
the diagonal matrix of eigenvalues of S and Q is orthogonal. Other representations may 
exist. For our purpose it is sufficient to pick any one of these candidates. The matrix A 
is only a tool with which to construct intermediate quantities and does not occur in the 
final transport equations, so this nonuniqueness does not lead to ambiguities. Given a 
choice of A, we may define standardized variables Zj(x), 



{{it- Xi){xj - Xj)) = Eij. 
In addition there may be nontrivial third moments, a^kit), 
((xi - Xi)(xj - Xj)(x k - X k )) = a ijk . 



(1) 



(2) 




(3) 
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The Zi have zero mean and orthonormal covariances. In this and subsequent expressions, 
the summation convention is applied to repeated indices. 

If the joint probability distribution of x is close to Gaussian (an assumption which 
we relax in Section [3]) , it may be represented by a Gauss-Hermite expansion 

a- 



P(x) d N x = P g {z) 



jk Hij k (z) 



d Jv x, (4) 



6 

where P g is a normalized Gaussian kernel, 

P « (z)= (det 2 U 1/2 eXP HO" (5) 
Note that Eq. (jlj) should be considered a function of x, via Eq. ([3]), although for 
convenience its right-hand side has been written in terms of z. Also, its nongaussian 
part has been written using a set of third-order moments af jfc associated with the Zj. 
These satisfy (ziZjZ k ) = a^ k and are related to the by the rule 

otijk = AiiAj m A kn af mn . (6) 

Thus, the a transform tensorially under the change of basis represented by Aij. The 
basis functions which underlie the cumulant expansion are products of Hermite 
polynomials. The n th polynomial in the Hermite sequence is obtained from Rodrigues' 
formula, 

H n (w) = (-l)"e- 2 / 2 £_ e -^/ 2 . (7) 
The H n satisfy an orthogonality relation, 

oo e -w 2 /2 



H n (w)H m (w) dw = n\5 mn . (8) 
The Hij k are defined by a generalized version of Rodrigues' formula. They satisfy 



H im .,„ = (-l)"«p -e W { -^) . (9) 

Eqs. (JHJ) (El) make the orthogonal in the measure exp(— z 2 /2) d^z. 

Eq. (jlj) is a cumulant expansion in the sense discussed in ^TJ It is numerically close 
to a Gaussian for small oc^ k . It has nonzero third moments, which are perturbatively 
small, but all higher cumulants are zero. Higher n th cumulants may be incorporated 
by including them in Eq. fl4]) as coefficients of n th order functions ifj r ..j n . Nevertheless, 
Eq. (jlj) is not an asymptotic expansion in the Z\. Moreover, finite truncations may be 
negative for some values of z,. For our purposes the cumulant expansion is a formal 
tool, and these mild pathologies are not a cause of serious difficulty because we do not 
make use of the probability distribution directly. 

The H iv .. in obey certain important identities. By repeatedly commuting z and d/dz 
one can use the generalized form of Rodrigues' identity to show 
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Further, differentiating Eq. and making use of fllOp one can show 

~p. = 3iimHi 2 ...i n + • • • + Si n7n Hi 1 ...i n _ 1 . (11) 

Eqs. ( TTOT) and (TTTj) play a significant role in extracting a moment hierarchy from the 
Gauss-Hermite expansion. 

2.2. Transport of the probability density 

Eq. (J3J) is time dependent, because of the explicit time dependence of X(£), and 
Oiijkif). We assume that time evolution of x is generated by a velocity field u(£, x), using 
the rule x = u. The vector u depends on x, but may also depend explicitly on time. It 
is possible to interpret u as a time-dependent Hamiltonian vector field whose integral 
curves are the allowed trajectories in phase space. As time evolves, the shape of the 
probability distribution is focused and sheared by the action of the velocity field. This 
geometrical evolution is described by the transport (or "continuity") equation 

dP diuiP) , , 

Eq. (Tl2|) accounts for changes in shape and profile of P, but conserves the overall volume 
of the distribution. It is the zero- diffusion limit of a Chapman-Kolmogorov or Fokker- 
Planck equation. These equations occur in many areas in physics, including the heat 
equation and Schrodinger's equation. 

In principle Eq. ( Tl2|) could be solved directly, but analytic progress is possible only 
for a limited choice of «j. Numerical approaches are rather involved. As an alternative 
to solving for P itself, Eq. (112j) may be converted into a system of coupled equations for 
the 77 th moments of P. The evolution equation for any moment will generally depend 
on all the others, but if increasingly high-order moments decrease in magnitude then it 
may be a reasonable approximation to truncate the coupled system at finite order. In 
what follows we will carry out this programme, assuming that it is only necessary to 
retain the moments Xi(t), £jj(t) and a^^t). 

In this section we are assuming that the third- and higher n th -order moments 
are perturbatively small. Therefore, probability is concentrated in the vicinity of 
the instantaneous centroid X{(t). The influence of the velocity field in reshaping the 
distribution is greatest in this region. Near the centroid, we find 

Ui = Uio + Uij{xj - Xj) + ^u ijk (xj - Xj)(x k - X k ) H 

= u io + u ijAj k Zk + -^Uij k Aj m A kn z m z n + • • • , (13) 
where the coefficients u^, Uij and Uij k satisfy 

■ (14) 



duAt) 

Uioit) = Ui(t)| x =X(t), Uij(t) - 



d 2 Ui(t) 

and Uijk{t) - 



=x(t) 



=X(t) OXjOXk 

A particle located at the centroid, X{ = Xi(t), would evolve according to dxi/dt = u i0 (t). 
The presence of the higher terms in (1131) reflects the way in which the wings of the 
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probability distribution, which are absent for a point particle, sample nearby parts of 
the velocity field. 

Whether or not depends explicitly on time, these coefficients will do so because 
they are evaluated at the time-dependent location x = X(t). Moreover, since we 
must truncate Eq. ( fTBT) at finite order to obtain a finite system of evolution equations, 
there will be an unaccounted remainder which makes the effective velocity field time 
dependent. 

To proceed, we must determine dP/dt and d(uiP)/dxi. According to Eq. ( Tl2l) their 
sum must be equal to zero. Using Eqs. ({TO]) and (TTTT) to exchange factors of z or d/dz 
(as applied to the H^...^) with sums of other i/- functions, it can be cast as a Hermite 
tableau of the form 

P g [cq + QHi + CijHij + c ijk H ijk H ] = 0. (15) 

The orthogonality of Hi v ..i n in the measure P g d N z implies a hierarchy of equations, 

c = Ci = c (i j) = c (i j fc ) = ■ ■ • = 0, (16) 

where brackets denote symmetrization of indices. 

The Co equation expresses overall conservation of probability, and is identically 
satisfied because Eq. ( fl2|) is conservative. The equation q = gives an evolution 
equation for the centroid, 

dXi 1 

-Q f -u i0 + -u imn Z mn + ---, (17) 

where '• • ■' denotes omitted terms which are higher order in cumulant expansion. The 
equation cnj\ = gives a similar evolution equation for the covariance matrix, E^-, 

d ^ij 1 1 

Finally, Cfe-H = is equivalent to an evolution equation for the 

= Ui m oi m jk + itj mri Sj m Sfc n + [i > j y k) + • • • , (19) 

where (i — > j — > k) denotes the preceding term with cyclic permutations of the indices. 
Eqs. (]17p -(ll9 p represent the first principal result of this paper. They agree with the 
corresponding evolution equations for first, second and third moments obtained in 
Ref. p], but apply for an arbitrary number of scalar fields. 



3. Method B: Generating functions 

In this section we are going to rederive the results of $2] using a technique which is valid 
for a probability distribution in the neighbourhood of an arbitrary kernel distribution. It 
can accommodate an arbitrary number of fields, and applies to any order in the cumulant 
expansion. The method applied in Ref. pQ and $2] exploited orthogonality properties of 
Hermite polynomials and was therefore restricted to a Gaussian kernel. Here, we take 
a different approach. We introduce generating functions for the cumulants, and show 
that this method leads to an efficient derivation of the evolution equations. 
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We retain the notation of §2J In this section, we introduce the new method by 
using it to study a probability distribution with a single field. The extension to an 
arbitrary number of fields involves more complicated combinatorics, but ultimately 
yields expressions equivalent to Eqs. (|T7|) -f fT9|) . We discuss the multiple-field case in 



Appendix A 



3.1. Generating functions 

In the one-field case, we denote the single field by x, and assume that the distribution 
of x is described by a time- dependent probability distribution of the form P(x, t)dx. 
We denote the mean value of x by X(t), so that 

*( t ) = /,P(,, t )d,. (20, 

In §2] we restricted our attention to the third moment, a. In this section we would like 
to generalize the analysis to all orders in the moment expansion. For this reason we 
introduce moments (i n (t), n = 0, 1, 2, . . . defined by 

fi n (t) = J [x-X(t)} n P(x,t) dx. (21) 

Since P is properly normalized we must have /io = 1, independent of time. Our definition 
of X implies fi\ = 0. Therefore the first nontrivial moment is the second, /i 2 - The entire 
infinite set of moments can be encoded using the "moment generating function" M(z, t), 
defined by 

M{z, t)= ! e< x -Vp(x, t)dx = jr Z -^M. (22) 

J n=0 U ' 

This definition ensures that the n th moment can be recovered using the identity 
8?M(0,t) =fin(t)- 

In $2] and the foregoing discussion we have restricted our attention to the moments 
of P. When discussing high orders in the moment expansion, it is sometimes more 
convenient to make use of the cumulants of P directly. We define the sequence of 
cumulants, their generating function C(z,t). This satisfies 

C(z,t)=\RM(z,t) = f^^&. (23) 

n=0 

Conversely, the moments are related to the cumulants through M(z, t) = e c ^ z,t \ Hence, 
knowledge of the moments is enough to determine all the cumulants, and vice-versa. 
For example, we always have = Ki(t) = 0, yU 2 (t) = ^(t), and /^(t) = ^(t), while 

m(t) = « 4 (t) + 3«|(t), (24) 
/2 5 (t) = K 5 (t) + 10« 2 (t)« 3 (t), (25) 
IMi(t) = «e(t) + 15« 2 (t)K 4 (t) + 10k§(<) + 15«l(t) (26) 

and so on ad infinitum. To obtain the cumulant K n one need only know the moments 
up to order n, and vice versa. 
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For n > 3 it is the n n which are most useful for the characterization of primordial 
nongaussianity. In the language of field theory, the cumulants are the connected 
correlation functions of x, while the moments are the disconnected correlation functions. 
In particular, a pure Gaussian distribution has only one nonzero cumulant, k 2 . On the 
other hand, all of its even moments /io, fJ>2, f>4, ■ ■ ■ , are nonzero. Thus, beyond third 
order, the cumulants provide a more suitable measure of departures from Gaussianity. 



3.2. Transport equation 

The transport equation is the one- dimensional version of Eq. ( fl2l) . Combining the 
transport equation with (|20|) yields 

dX f d f 

= x —P(x,t) dx = u(x)P(x,t) dx (27) 
at J at J 

As in £j2j we assume that x evolves under the influence of a velocity field which can be 

expanded around the instantaneous centroid X(t), 

J2^1[ x -X(tT, (28) 



u[x 

n=0 



where, in comparison with Eq. (TIB"]) , we have adopted a slightly different notation in 
which u n denotes the n th derivative d n u/dx n . Returning to (127]) and applying (128]) 
together with the definition of the moments, we find 

ix = ^_jf^ = uo(t) i (()M2(() + ... (29) 

n=0 

With the evolution of X(t) in hand, we can derive the evolution of all other 
cumulants. Using ( 1231) . we conclude that the time derivatives dn n /dt obey 

Ez n dn n {t) _L y-v z n dfx n (t) 
n\ dt ~M(z,t)^n\ dt ' 1 ' 

n=0 v ' ' n=0 

from which it follows that dC/dt is the generating function for the cumulant time 
derivatives. Following steps similar to those which led us to fT27|) . it can be seen that 
the evolution equation for each moment can be written 

dn n (t) 



dt 



^ -Jj ^ji n +k-l(t) ~ fl n ^l(t)ll k (t) U k (t) (31) 



k=0 

Inserting this expression in ( 130]) yields a generating function for the time derivatives of 
each cumulant. The result can be written in terms of the moments. Finally, using the 
relationship between the moments and cumulants, the cumulant evolution equations can 
be written in terms of the cumulants alone. This gives a closed system of equations for 
the cumulants, with an infinite number of variables. 

In practice, only a finite number of variables can be evolved and it is necessary to 
truncate both the series of cumulants and the expansion of u. If the Taylor expansion of 
the velocity field is truncated, then it is evident from ( 13~TT) that the evolution equation 
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for each moment involves moments of higher order. Accordingly, the time-evolution of 
an individual cumulant always involves cumulants of higher order. Thus, even when the 
velocity field is truncated, the system of cumulant evolution equations does not close at 
any finite order. We must choose an order at which to approximate the full cumulant 
expansion. If desired, we can obtain an estimate of the error involved in truncation at 
this order by determining the degree to which time evolution sources higher cumulants, 
forcing them to become nonzero. 

We now illustrate this procedure in operation using a simple example. Suppose we 
wish to carry the velocity field expansion to third order in (x — X), so that 



u{x) = uo(t) + [x- X(i)]ui(i) + y i x 



X(t)Yu 2 (t) + -[x 



Assume further that the maximal nonzero cumulant is of order four 
formulae above, we find 

= «o(*) + ^U 2 {t)K 2 {t) + ^M 3 (*)« 3 (*) + ^ [^ 2 (tf 

The cumulant evolution equations enforce dno/dt 
the cumulants k 2 , k 3 and k 4 evolve according to 



X(t)} 3 u 3 (t). (32) 
Applying the 



(33) 



dn 2 
~dt 

dt 
dft 4 

~df 



2UxK 2 + U 2 K 3 + U 3 

3mik 3 + u 2 

4uiK4 + U 2 



K 



1 

2 + 



3 Ate 



2* 



-jU 3 K 2 K 3 



12k 2 k 3 + 2k 5 



u 3 



and dfti/dt = 0. Furthermore, 

(34) 
(35) 
(36) 



4^2 + 6K3 + 



where we have suppressed the time dependence of the Kj and uj. 
Under our assumptions, the evolution equation for k 5 is 
d« 5 



dt 



It 2 



20^2^4 



154 



u 3 



30^2^3 + 25k 3 k 4 



(37) 



If a truncation to fourth-order quantities were self-consistent, this equation should vanish 
as a consequence of our assumption that n n = for n ^ 5. That it does not vanish 
is a measure of the error involved in our truncation. As explained above, a similar 
effect occurs no matter at which order the truncation is made. In the analysis of this 
section, we have made no assumption that either the moments or cumulants are small, 
or order themselves into an ultimately decreasing sequence. Eq. (1371) demonstrates that 
the truncated hierarchy will be a useful predictive instrument only when a negligible 
variation in k 5 is sourced over the time interval of interest. This will typically (but 
not absolutely) require the cumulants to fall in an ordered structure, for example 
\k 2 \ > \k 3 \ > \k&\. 

For example, we may suppose that \n n \ = 0(S n ) where 5 1 is a small positive 
number. Eq. (1371) shows that source term for k§ is of order 0(5 6 ), and therefore after a 
short time interval At we can expect k 5 ~ (At)<5 6 . Truncation to fourth order, setting 
«5 and all higher cumulants to zero, may be acceptable approximation over sufficiently 
short times that k 5 does not grow to the degree that it contaminates any observable 



Moment transport equations for the primordial curvature perturbation 



11 



of interest. How long this time can be is model dependent. A similar analysis can 
be given for all higher cumulants. Growing secular terms which eventually invalidate 
perturbation theory are a typical feature in studying the evolution of fluctuations. They 
are encountered directly in the Feynman amplitudes of quantum field theory, and in the 
most common implementations of the separate universe picture. 

Eqs. f l34"|) - fl36l) coincide with Eqs. ffTB"]) and ffT9l) when specialized to the single- 
field case, but include the contribution of the third derivative of the velocity field, u 3 . 
Eq. (1361) gives, for the first time, the evolution of the kurtosis in a single-field setting. 
In Appendix A we apply this method to the case of multiple fields, and again find 
agreement with the results of §2.11 



4. Inflationary perturbations 

We now wish to apply the general framework assembled in §§2H3] to the fluctuations 
which are generated and subsequently evolve during an inflationary era. To do this we 
identify the variables Xi of §§2H3] with the light, scalar degrees of freedom which are 
excited at that time. We must also identify a choice of time variable, labelled t in §§2rEI 
which corresponds to a choice of slicing in the spacetime picture. 



4-1- Gauge choices 

Consider a theory of M scalar fields coupled to a metric theory of gravity. The degrees 
of freedom in this system are the M fields themselves, <fti for i e {1, . . . , M}, together 
with the lapse and shift functions, written N and N a , and the spatial 3-metric h a b, 

ds 2 = -N 2 dt 2 + h ab {dx a + N a dt){dx b + N b dt). (38) 

In Einstein gravity the lapse and shift are not propagating fields, and are eliminated by 
constraint equations. The 3-metric h a b encodes six independent degrees of freedom but 
in this paper we concentrate on the two spin-0 modes, of which only one is physical. 
We take this to be the volume modulus det h a b\ the other can be absorbed by a spatial 
coordinate redefinition. We write det h a b = & 6 ( N ~ N °'8 a b, and refer to N as the integrated 
e-foldings of expansion. The zero point N is arbitrary. Note that the integrated 
number of e-folds is quite separate from the lapse function which occurs in Eq. (|38|) . 
also traditionally denoted N. In what follows we shall work in an unperturbed universe 
for which the lapse is unity. Therefore we have no need to refer to it explicitly, so that 
the appearance of N without qualification is unambiguous. 

We now have M + 1 variables: the M fields, and the integrated expansion N. Not 
all these are independent, and one of them can be eliminated by a suitable choice of 
time, t. Whatever choice we make, surfaces of constant t foliate spacetime into spatial 
hypersurfaces which we refer to as a slicing. The transport equation, Eq. ffl2l) . evolves 
a probability distribution for the Xi from one slice in this foliation to the next. There 
are two slicings of particular importance for inflationary fluctuations. 
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Spatially flat slicing. In this gauge, spatial hypersurfaces are chosen so that the 
integrated number of e-foldings, N, is uniform across the slice. The scalar fields fa 
fluctuate from place to place. To apply the formalism of £[2]we work on slices of uniform 
expansion N, and smooth the scalar fields on a comoving lengthscale L to give an 
ensemble of L-sized regions. 

These regions traverse a bundle of adjacent trajectories in field space, with some 
characteristic dispersion and higher-order moments which it is our intention to calculate. 
To make contact with observation the bundle should be chosen so that every trajectory 
reheats almost surely in our local vacuum. With this choice we can suppose the 
dispersion and higher-order moments of L-sized regions in our observable universe should 
be similar to those of the bundle. On scales comparable to the size of our observable 
patch we sample only a small number of trajectories, making the prospect of a mismatch 
with the bundle average ("cosmic variance") more likely. 

The centroid of the bundle follows a path $f(t) in field space, where the superscript 
'L' indicates dependence on the smoothing scaleqj The probability distribution we wish 
to calculate is a function of the scalar fields, fa = $f (t) + 5<pf. We take Xi = <pi/Mp and 
Xi(t) = $f (t) /Mp. The probability density obtained in this way gives information about 
the distribution of field values on the scale L only. To obtain a relation between bundles 
smoothed on different scales L and L', two separate distributions must be computed 
and their information combined. This would be necessary, for example, to obtain the 
spectral index. 

We also require initial conditions, set at the time ti when the wavenumber 
corresponding to the smoothing scale L crosses the horizon. We pick initial expectation 
values $f (tz,) centred on the inflationary trajectory of interest. In the spatially flat 
slicing, the joint probability distribution of fluctuations in the scalar fields at time t^ 
can be calculated directly [9], [101 HH [12] , with each field acquiring a variance of 
order (5<j)f /Mp) l = 0"£, where o~l = Hi/Mp ~ 10~ 5 and Hl is the Hubble rate when 
the wavenumber 1/L crosses the horizon. The skewness of the bundle is negligible, 
with a^ k ~ a\ ~ p31 [2]. In addition, to leading order in e = —H/H 2 , the slow- 
roll condition makes the fields uncorrelated at horizon exit, with (5<fri5<j)j/Mp)L ~ e L0~\ 
if i 7^ jlj] Taking the initial dispersion to be given by ctl and setting afj k and any 
cross-correlations to be initially zero, we can evolve the a and a from horizon crossing 
until any desired future time, such as the end of inflation. At this point the field's 
moments on the flat hypersurface can be used to calculate the moments of ( using a 
gauge transformation. For two fields this proceedure was performed in Ref. pQ, and the 
technology developed in this paper makes this possible for an arbitrary number of fields. 

f Note that &f (t) need not be an integral curve of the velocity field, and therefore may not constitute 
an allowed inflationary trajectory. 

| The conclusion that the 8<j)i are virtually uncorrelated at horizon exit implies that the inflationary 
trajectories — described by the integral curves of the velocity field in this gauge — are effectively 
straight lines for a few e-folds around horizon-crossing, up to corrections of 0(e). This applies in 
canonical inflationary models, but in more general examples this may not occur [14j . In such cases the 
theory becomes more complicated and the formalism of this paper will no longer apply. 
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Uniform density slicing. Alternatively, we may choose our spatial hypersurfaces so 
that the Hubble rate, H, is uniform over the slice. In Einstein gravity the Friedmann 
constraint enforces 3MpH 2 = p, so slices of uniform H are also slices of uniform 
density. On this slicing the integrated number of e-foldings, N, will typically vary 
from place to place. The counting of independent fluctuations is the same as in the 
spatially flat slicing, because the Friedmann constraint makes one field a function of 
all the others and the Hubble rate, H. Without loss of generality we can suppose that 
(pM = 0m(0i, • • • , 0a/-i, H). For notational convenience we define Sj = fa/Mp. To apply 
the framework of §£j2]-[3]we must set the time variable, t, to equal H/Mp and choose the 
variables Xi whose distribution we wish to calculate to be {s\, . . . , sm-i, N}. 

It is still necessary to choose a smoothing scale, L. The centroid of the 
bundle is characterized by the expectation values of the first M — 1 scalar fields, 
{3>i(£), • ■ • , ^M-i(^)} together with the mean integrated expansion experienced by 
trajectories within the bundle, Ni,(i). We write N = N L (t) + ( L , where Cl is the 
uniform density gauge curvature perturbation smoothed on scale L. On the large scales 
we are considering, the uniform density gauge and comoving gauge coincide which makes 
(l numerically equal to the comoving curvature perturbation, TZl- 

4-2. Moment transport in the uniform density slicing 

To implement moment transport in the uniform density slicing, we will require initial 
conditions for the dispersion and higher moments of the fluctuations {s\, . . . , Sm-i, C}- 
These have not yet been calculated directly, but can be obtained from the joint 
probability distribution in the spatially flat slicing [2] in conjunction with gauge 
transformations relating super-horizon quantities in the uniform density gauge to those 
in spatially flat gauge. These transformations can be obtained using coventional 
cosmological perturbation theory [15] , or using the separate universe assumption, which 
was the approach taken in Ref. pQ. Using the separate universe assumption, one writes 
C using the '<5iV formula [16] 

CO) = N(t, fa + 5 fa) - N{t, fa) = NiSfa* + -NijSfaJfa* + ■■■, (39) 

where N(t, 0*) measures the e-foldings between a spatially flat slice on which the fields 
take prescribed values fa* and a subsequent uniform density slice at time t. We are 
considering the special case of flat and uniform density slices which coincide on average, 
and are therefore only perturbatively separated. The coefficients satisfy iVj = dN/dcf)^ 
with similar definitions for iV^ and higher derivatives. 

Initial conditions for isocurvature fields. The fluctuations {si, . . . , sm-i} can be 

calculated using an analogous formula, as described in Ref. [I], 

dfa 1 d 2 fa 

Si Mp = T^Sfa* + - I 1 5fajfa* + ■ ■ • , (40) 

The superscript 'c' denotes scalar fields evaluated on a comoving (uniform density) 
spatial slice, in the same way that V denotes fields evaluated on spatially flat slices. 
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Eq. gives each field Sj a dispersion erf, approximately satisfying 



E^T^ (no sum on z) (41) 



1- M| — „..,,. 



together with negligible third moments and cross-correlations. By a suitable choice of 
origin we can arrange that N = on the initial spatially flat slice, and therefore N 
satisfies 

on the initial uniform density slice. Since this is very small it is a reasonable 
approximation to take Nl ~ 0. Eq. ( |39|) generates nonzero correlations between £ 
and the Sj, allowing us to calculate the covariances (C s «) l- It can also be used to obtain 
the ^-dispersion, (£Q L . We find 

Expressions for velocity field. With this choice of variables and the assumption 
that third-order moments are zero initiall}{§|, the initial conditions, Eqs. ( TT7|) - (|T9|) . can 
be used to compute the covariance matrix and third-order moments at any later time, 
provided expressions can be found for the velocity potential and its derivatives. In the 
uniform density gauge, the allowed trajectories are integral curves of 

, , dsj , , , dN . 

u i = M P^TT and u n = M p — . 45 
an aH 

Applying the slow-roll approximation allows us to write the velocity potential as a 

function of s, and H, and we find 



and 



u, = -f ^ (46) 

In order to write these formulae we have defined a set of partial slow-roll parameters, 
6j, which satisfy 

. _ 1 # 



2MIH 2 ' 



(4* 



§ It is clearly possible to extend the approach we have described to calculate the initial third order 
moments in the uniform density gauge, given the known initial conditions in the flat gauge. This 
requires the use of N^j and d 2 4>i/d<f>jk*- The assumption that these are zero, however, introduces an 
error only of the same order as that already present from the typical approximation of taking the initial 
third moments to be zero in the flat gauge. 
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where <fi = dV/d<pi/3H . To leading order in the slow-roll approximation, e = £\ Note 
that both Ui and are independent of N, but depend explicitly on the time variable 
H . Therefore, all ^-derivatives of the velocity field vanish identically. Derivatives 
of Eqs. (|46|) and (|47|) with respect to the Sj can be found after using the Friedmann 
constraint to eliminate the variations d<j)M/d4>i and d<pM/dH. 

4-3. Numerical performance of the moment transport method 

The imminent arrival of high-quality microwave background data from the Planck 
satellite implies that it will soon be necessary to obtain accurate estimates of the 
nonlinearity parameters in a wide range of inflationary models. Although analytic 
formulas for /nl exist, they are available only for certain forms of the potential and even 
when their use is possible they rapidly become unwieldy in the limit of a large number 
of fields. For these reasons we expect numerical methods for computing /nl to become 
of increasing importance. In §5] below, we discuss a numerical implementation of the 
moment transport method using the uniform density slicing. Here, we briefly comment 
on the computational efficiency of the moment transport algorithm in comparison with 
alternative approaches. 

In a system with M fields, the moment transport method requires a solution of 
Eqs. (|T71)-(|T9l). These comprise M unique equations for the centroid, M{M + l)/2 for 
the covariance matrix and M(M + 1)(M + 2)/6 for the 3-point functions. Therefore, 
for large M we must solve a coupled system of 0(M 3 ) equations. In comparison with 
popular alternative methods based on the 5N formula [T71 [181 CHI EES], we argue that 
the transport method is computationally simpler. 

To make use of the SN formula requires calculation of the derivatives dN/dfa*, 
d 2 N/d<f)i*d(f)j*, and so on. Therefore a direct implementation of the 5N formula 
with M fields requires numerical evolution of the background field equations for many 
initial conditions, from which the necessary derivatives may be extracted. 

How many evolutions of the background equations are required? This will determine 
the computational efficiency. To calculate /nl, we require derivatives up to second order 
but not higher. The first-order derivatives may be obtained by taking finite differences 
between inflationary trajectories with initial conditions separated by a small distance 5, 
giving derivatives up to an error of order 5 2 . This requires of order M + 1 evolutions of 
the M equations, or the solution to ~ M 2 ordinary differential equations!]]] Extending 
this argument to the second order derivatives shows that to determine /nl we must solve 
~ M 3 ordinary differential equations. For higher moments the same counting applies, 
so that an evaluation of the trispectrum will typically require the solution to ~ M 4 
ordinary differential equations, or more generally ~ M n for the amplitude of n-point 
correlations. 

We conclude that, to compute /nl, both the transport method and a direct 5N 
require the solution of 0(M 3 ) ordinary differential equations. Asymptotically, their 

|| In practice, a more accurate discretization scheme may be required. 
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relative efficiency depends on details of the algorithm. However, the method of moment 
transport is especially simple. Numerical 5N requires a discretization scheme to compute 
the partial derivatives. The final accuracy can depend on our choice of discretization. 
Also, as we will discuss below, SN requires the e-folds of expansion to be determined very 
accurately on successive time slices. In comparison, the transport method requires only 
the solution to a set of coupled ordinary differential equations. This allows off-the-shelf 
differential equation solvers to be brought to bear on the problem immediately. 

In a direct SN algorithm, very high accuracy is required when evolving the 
background field equations because at the end of the calculation we must take finite 
differences to construct the derivatives of N. In simulations with up to M = 5 fields 
we have found this process to be sensitive to small numerical inaccuracies; the moment 
transport algorithm produces accurate results with larger numerical tolerances. This 
observation can be explained in a simple way. Suppose we wish to compute SN to 
fractional accuracy /, where 

/ = (49) 

On the one hand, if we use the naive SN algorithm and compute N using an integration 
routine which operates at fractional accuracy fx, then the absolute error in N will 
be fxN. Therefore, the absolute error in 5N is also f\N, and to achieve the target 
precision (I49p we must choose f\ ~ f5N/N. On the other hand, the moment transport 
approach essentially integrates SN directly. Using an integration routine with fractional 
accuracy f 2 we will evaluate 5N with the same fractional accuracy, so / 2 ~ /. Since 
5N/N < 1, we have fx <§C f%- We conclude that the moment transport algorithm can 
operate with much lower numerical tolerances than the naive 5N approach. 

Finally, we note that direct implementation of 5N is not the only possibility. 
Yokoyama, Suyama & Tanaka [201 EE) suggested an approach which is broadly similar 
to that employed in Ref. [I], reviewed in Ref. [22]. In this approach, one seeks to 
calculate the field perturbations on a uniform curvature hypersurface at the time of 
interest, as functions of their initial values, using the 5N formula to effect the final 
gauge transformation. One important difference between the formulation of this paper 
(and Ref. PQ) and that of Yokoyama et al. is that the authors of Refs. (20], [2TJ continue 
to work in terms of the field perturbations, rather than evolving the moments of the 
distribution. We believe this makes our formulation simpler to implement in practice. 



5. Numerical examples 

In this section we present results obtained from the moment transport method, 
implemented on a desktop computer. This demonstrates that it is possible to compute 
/nl m systems with at least M ~ 10 2 fields using commodity hardware. The exact time 
required for a simulation is highly model dependent, but a simple M — 10 example may 
only take a few seconds, while a M = 100 calculation can be carried out in less than 
an hour. Our numerical code is implemented using Matlab, and evolves the uniform 
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density gauge transport equations, Eqs. (TTTT ) — (TTOT) . Using this code we are able to obtain 
numerical solutions for the potential [231 12H I2S1 [26] 



This describes a number of uncoupled fields with quadratic potentials. It is also the 
small-field approximation to a collection of uncoupled axions, which have trigonometric 
potentials; this latter case is often known as Nflation. The variables evolved in this 
gauge and their initial conditions are determined by applying the discussion of §4.21 
In the special case where all are equal there is an O(M) symmetry, and the fields 
roll radially to the origin. Where this symmetry exists we have verified that our code 
reproduces the expected single- field result of constant ( and negligible /nl, for up to 
10 2 fields. 

Where the m; are different, it is known from analytic calculation that this potential 
does not give rise to a large nongaussianity [27] [28]. Our results confirm this conclusion, 
but the simplicity of the model and the existence of analytic predictions makes it a 
useful test of our method. In Figs. Q~H3]we show illustrative results. Consider first the 
case of a small number of fields. Rigopoulos, Shellard & van Tent [291 EH] analysed a 
two-field example with mass ratio mi = 9m>2. Vernizzi & Wands later studied the same 
model using a combination of analytic and numerical methods [27] and concluded that 
/nl at the end of inflation was very small, of order 10~ 2 . 

We depict the evolution of /nl in Fig. [TJ For comparison, we plot a calculation 
of /nl, for the same model and initial conditions, using a numerical slow-roll 
implementation of the SN formula, finding complete agreement. In this and subsequent 
figures the horizontal axis is labelled by N, the number of e-folds which have elapsed 
since horizon-crossing of the wavenumber of interest. Because our moment-transport 
and SN calculation are performed using the slow-roll equations of motion, we do not 
need to allow a 'lead time' for the simulation to converge to the slow-roll attractor. 

The qualitative behaviour of /nl i n this model was discussed at the end of §4.21 
Our evolution agrees with this discussion, and also the explicit calculations of Vernizzi 
& Wands [27] and the moment transport method using the spatially flat slicing [1] . The 
most prominent feature is a well-documented spike, which occurs when the heavier field 
reaches the vicinity of its minimum and decouples from the dynamics. 

In Fig. [21 we increase the number of fields to 10. For simplicity we take fa = 5 
for all fields and distribute the masses logarithmically, in such a way that m; = 2mj_i. 
As successive fields evolve to their minima, a wavelike structure is produced in /nl- Its 
value at the end of inflation is again of order 10~ 2 . Kim & Liddle [26] argued that, at 
the end of inflation in a model of the form ( )50|) . the bispectrum imprinted for a mode 
of _ber k could be parametrized * » 1/JW.fl In this estimate, N m 

f We have neglected a quantum-mechanical contribution generated by interference among field modes 
at horizon crossing. This contribution is not determined by the moment transport method, and its 
contribution is not represented in Figs. fTHU It is known to be small [31] [27] and may be neglected 
when |./ NL | > 1. 




(50) 
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Figure 1. Evolution of /nl (solid red line) calculated using the moment transport 
equations for a two-field Nflation model, working in the quadratic approximation. The 
green dashed line shows /nl calculated for the same model using the 5N formula. 
Initial conditions are described in the text. In this and subsequent figures the horizon 
axis is labelled by N , the number of e- folds which have elapsed since horizon-crossing 
of the wavenumber of interest. 
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Figure 2. Evolution of /nl calculated using the moment transport equations for a 
10-ficld Nflation model in the quadratic approximation. The initial conditions and 
distribution of masses are described in the main text. 




Moment transport equations for the primordial curvature perturbation 



19 




N 

Figure 3. Evolution of /nl calculated using the moment transport equations for a 
100-ficld Nflation model in the quadratic approximation. The initial conditions and 
distribution of masses are described in the main text. 

is the number of e-folds to the end of inflation from the field values at horizon exit of 
mode k, dropping any correction from the end of inflation. Within this approximation, 
the result is independent of the masses, initial conditions, or number of fields. In this 
model, we find iV* « 63 yielding a value for /nl in good agreement with the formula of 
Kim & Liddle. 

In order to demonstrate the ability of our algorithm to deal with a large number 
of fields, we give two examples with M = 10 2 fields. First, we retain the potential (150]) 
and set <pi = 1.5Mp for each field. We distribute masses such that nii = mj_i + O.lmi. 
The evolution of /nl in this model is shown in Fig. [3] At the end of inflation, /nl is 
marginally smaller than in the 10-field case. This suppression can be thought of as a 
consequence of the central limit theorem, which requires that ( becomes Gaussian if it 
receives comparable contributions from a large number of field perturbations with finite 
variance. To achieve a large nongaussianity as M becomes large, one must arrange that 
( be dominated by only a few fields [32], as we will discuss below. 

These plots show that the horizon-crossing approximation is valid, in these models, 
to within roughly 5% and 15%, in the 10- and 100-field cases respectively. In the latter 
case, this discrepancy can likely be ascribed to the relatively large fraction of fields 
still in motion at the end of inflation. As described above, the typical effect causes 
/nl to decrease as a consequence of the central limit theorem. As the horizon-crossing 
approximation begins to fail the value of /nl ceases to be universal and acquires a 
dependence on the details of the model. 
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Second, we study /W i n an Nflation model which retains the full trigonometric 
form of each potential, |jj 



It has recently been shown that these models have a phenomenology quite different from 
Eq. (150]) if the initial conditions populate regions in field space close to the maximum of 
the potential, where the quadratic approximation is poor [32]. We choose the fi to have 
a common value, /. Where N fields explore the hilltop region and contribute roughly 
equally to the curvature perturbation, Kim et al. [32] estimate that it is possible to 
achieve /nl of order 



This can easily be of order 1-10 for / « Mp and only a single field in the vicinity 
of the hilltop. Therefore, Eq. (I5ip constitutes an excellent test that our methods are 
effective in models for which /nl is not negligible. 

We work with M = 10 2 fields as before, and set / = 5Mp. We fix initial conditions 
so that <pi = 1.25Mp for all but one field, and assume this remaining field is very close 
to its potential maximum, with <p = 2.49Mp. Under these conditions Eq. fl52|) yields 
/nl ~ 0.66. We plot the evolution of /nl in this model in Fig. HI Initially, the field 
closest to the hilltop is held in place by the large Hubble friction generated by all the 
other fields. This is the classical assisted inflation mechanism [33] . Our choice of fi and 
initial conditions implies an 0(M — 1) symmetry among the remaining fields, which roll 
radially away from the hilltop. In a more general model, the fields furthest from the 
hilltop would be sequentially ejected into their minima, where they decouple from the 
dynamics. During this phase /nl is constant and practically zero. 

Eventually, the Hubble friction decreases sufficiently to allow the field closest to 
the hilltop to roll. While this field is near the maximum of the potential it can support 
a few e-foldings of inflation, but it is rapidly ejected from the vicinity of the hilltop 
and accelerated expansion ceases. During this process, /nl suddenly receives a large 
contribution from the latent nongaussianity which was imprinted in the fluctuations of 
this field around the time of horizon crossing^ This contribution to /nl is proportional 
to the curvature, or //-parameter, of the cosine potential near its hilltop region. As the 
single hilltop field joins and then decouples from the dynamics, there is a transient spike 

| We carry out this calculation in the spatially flat slicing, which validates our M -field formulae in this 
gauge. 

§ This illustrates a subtle feature of the 'horizon-crossing' approximation, used in Refs. [26J [32] - If all 
trajectories converge to an attractor, then the statistics of the curvature perturbation are determined 
only by the fields' values at the time of horizon crossing. In the horizon-crossing approximation, only 
this contribution is kept. This does not imply that the final /nl has any relation to its value actually 
at the time of horizon crossing, and indeed it will typically be quite different. In the present model 
/nl ~ for almost the whole history of inflation, only achieving its 'horizon-crossing' value as the 
inflationary phase comes to an end. 




(51) 




(52) 
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Figure 4. Evolution of /nl calculated using the moment transport equations for a 
100-ficld Nflation model, retaining the full trigonometric form of the potential. The 
initial conditions are chosen so that a single field explores the region close to the hilltop. 
The nongaussianity of this field dominates the late-time attractor solution, at which 
/nl converges to a time-independent nonzero value. 



where /nl grows rapidly, and oscillates in sign. At the end of inflation, it settles down 
to a time-independent value /nl ~ 0.65, which is an accurate match to the prediction 
of Kim et al. [32], obtained using the horizon- crossing approximation. 



6. Conclusions 



In this paper we have refined and extended the "moment transport" method, introduced 
in Ref. P to compute the bispectrum nongaussianity parameter /nl in a two-field model 
of inflation. The formulation given in this paper contains three significant improvements. 

First, the results of Ref. [1] were valid only for a two-field model. We have 
extended the evolution equations quoted in that paper to an arbitrary number of fields, 
confirming the conjecture made there that these equations were field-space covariant. 
This is essential if our formalism is to be applied to models with a large field content, 
such as Nflation. It is also a practical requirement in many quasi-realistic models of 
inflation — perhaps deriving from supergravity, string compactifications, or the scalar 
fields available within the MSSM — which typically contain a more modest number of 
fields of order M ~ 10. 

Second, we have shown how to write the evolution equations for the curvature 
perturbation ( directly. As an ancillary benefit this leads to a simple formula for the 
time evolution of the /nl parameter. Third, discussed in §§3]and Appendix A , we have 
developed an entirely different method of deriving the moment transport hierarchy, 
making use of cumulant expansion in the distribution of field values. The original 
version of the moment transport method made the simplifying assumption that the field 
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distribution was nearly Gaussian, by expanding the true distribution as a perturbative 
series around a Gaussian "kernel" distribution. Our new technique works with the 
probability distribution directly and does not require a Gaussian kernel, providing an 
alternate derivation of the moment transport system and extending its range of validity 

Both versions of the moment transport method share the advantage that they are 
simpler to implement numerically, and may have accuracy and performance advantages 
over a direct 5N algorithm. Similar properties may be shared by the approach of 
Yokoyama et al. [201 ETJ- This is because they require only the solution of a set of coupled 
ordinary differential equations. Furthermore, error propagation in the moment transport 
system can be controlled more straightforwardly, since one evolves the nongaussian 
moments directly and need not compute differences of large quantities, as in the SN 
framework. The ( version of the moment transport procedure presented here also 
highlights the source terms which are responsible for the growth and decay of the 
moments of (, and of /nl- We hope that this will eventually clarify the origin of 
nongaussianity on a model-by-model basis. 

Our method is limited by the requirement that the cumulant expansion remains 
valid at all intermediate times. We may expect a breakdown when the third-order 
cumulants become comparable to the variance, so as a conservative estimate we expect 
the expansion to be valid provided /nl < 10 4 throughout the evolution. A similar 
limitation is shared by any approach which uses the evolution equations of perturbation 
theory [201 [2TJ , but is avoided by direct numerical SN. 

In our method, there is no obstruction to dropping the inflationary slow-roll 
assumption, except at horizon crossing where it is needed to fix the initial conditions for 
the moment transport equations. We have made use of the slow-roll approximation in 
our numerical computations in Ref. [1] and $51 Relaxing this assumption would entail 
additional initial conditions and evolution equations for the field velocities. This can be 
accomplished quite naturally in the moment transport framework, since the second-order 
differential equations describing the non-slow-roll dynamics can always be expressed as 
a larger set of first-order differential equations. This could be achieved by doubling the 
variables Xj of $2] and identifying the new variables as field velocities. The equations of 
motion would reduce to a velocity field on this doubled space, and our formalism would 
go through as before. 
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Appendix A. Multiple field evolution equations 

The derivation of the evolution equations for multiple fields follows essentially the same 
steps as the one-field case, but with more complicated combinatorics. We assume we 
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have D fields, with probability distribution 

P(xi, x 2 , ...id, t) dxidx 2 ...dxD = P(x, t) d D x 
We define the mean position of the distribution X, by 

and denote the moments by 

Vmn 2 ...n»(t) = j (x 1 -X 1 (t))^(x 2 -X 2 (t)) n \..(x D -X D (t)) n °P(x,t)d D x (A.3) 

The rank of a given moment ^mn 2 ...n D {t) is n\ + n 2 + ■ ■ • + n£>. The moment generating 
function M now has D dummy variables Zj, and is given by 
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The cumulants are defined by the cumulant generating function 
C(z u z 2 ,...z D ,t) = \nM(z 1 ,z 2 ,...z D ,t) 
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The cumulants of zeroth, first, second, and third rank are identical to their corresponding 
cumulants. At higher order, we have for example (when D = 2) 

- 2 (A.8) 

(A.9) 
(A.10) 
(All) 

yU 4 = «04 + 3Kq 2 (A. 12) 

The equation for the cumulant time derivatives is then analogous to the one-field case. 
For example, using the multi-field probability conservation equation 

^^ + tirl^)PM]=0 (A.13) 
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8xa 



we find 



dX, 
dt 



XiP(x) d D x 



u,P(x) d D x 



If we define the velocity field expansion in the natural way by 

00 

^ J j\nin2...ri£) 



E 



ni,ri2...ri£)=0 



ni\n 2 \...nr>\ 



x 1 -X 1 ) n Hx 2 -X 2 ) n \..(x D -X_ 



D 



(A.14) 



(A.15) 
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then we find 

dt ^ n l \n 2 \...n D \ 

m,n2...rt£)=0 

which is clearly the multi-field analogue of ( 1291) . 

Just as in the single-field case, we can use (1A.16|) to derive the equations of motion 



for the cumulants. As in the single-field case, we have 

„nl „n2 7 nfl J 

^2 ■■■ /C D u - h >n 1 n 2 ...n D _ /» „n 

niy...iin! dt 

ni,n,2...n D =0 

I ~nl „n2 ,nfl J ,, 



M(zi, Z2, ...ZD,t) ^ n x \n 2 \...n D \ dt 

Following a derivation which entirely parallels the single-field case, we arrive at the 
expression 

dfI ni n 2 ...n D n j(^-j + Bj) / a i q\ 

dt ^-^ / — ' m 1 !m 2 !...mn! 

mi,m2,...m£,=0 j=l 



where 



and 



P"ni+mi,n2+m2,...nj+rn.j — l,...ri£)+'mD (A. 20) 



Bj f^ni,n2,...nj—l,...n£)f^mi,m2,...m£) (^"21) 

These equations, combined with the expressions for the moments in terms of the 
cumulants, enables the system of equations for the cumulants to be derived. 

The algebra involved in deriving the evolution equations for the cumulants is 
straightforward, but tedious. Fortunately the required manipulations are entirely 
mechanical and can be implemented in a computer algebra system, such as Mathematica. 
As an example, if we expand to the third cumulant and to quadratic order in the velocity 
field, we find 



dn 



20 



dt 

d/cu 
dt 



dn 



02 



and 



dt 

d^30 

dt 
d/c 2 i 
dt 



2lil|oiKll + 2Wi| 10 K20 + Ml|02^12 + 2li 1 | 11 K 2 l + «1|20^30 (A. 22) 

Ml|01^02 + Ul|10«41 + U 2 \0lK U + M 2 |10^20 + Wl|ll«12 

11 11... 
+ —^1|02^03 + -^2|02^12 + U2|ll«21 + ^^O^l + -M 2 |20«30 (A. 23) 

2«2|01^02 + «2|10«11 + ^2|02«03 + 2u 2 |ll«12 + n 2 |20^21 (A. 24) 



3wi|02«ii + 6iti|iiKnK 2 o + 3mi| 20 ^2o + 3w i|oi«2i + 3mi|i k 30 (A. 25) 

2^1|02^02^11 + 2«i|ii« n + U2|02«n + 2ui|01«i2 + 2u 1 | 11 k 02 «; 2 o 
+ 2Mi| 20 KhK20 + 2u 2 |ilKuK 20 + M2|20^20 + 2 Wl|10«21 
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d«i2 
dt 



+ «2|01 K 21 + «2|10 K 30 (A. 26) 

U 1|02 K 02 + n l|01 K 03 + 2u 1 | 11 Ko2 K ll + 2n 2 |02^02^11 + n l|20 K ll 



+ 2u 2 |ll/«n + Wi|i Kl2 + 2m 2 |01^12 + 2« 2 |iiK 2/S20 

+ 2m 2 | 20 kiik 20 + 2w 2 |i k 2 i (A. 27) 

— ^- = 3m 2 | 02 k 02 + 3m 2 | iKo3 + 6m 2 |hk 02 kii + 3tt 2 |20«ii + 3w 2 |i ki 2 (A. 28) 

These equations can be generalized to any number of fields, or any order in the cumulant 
expansion, by using the expressions given above. 
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